Iterative summation of path integrals for nonequilibrium molecular quantum transport 
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We formulate and apply a nonperturbative numerical approach to the nonequilibrium current, 
I{y), through a voltage-biased molecular conductor. We focus on a single electronic level coupled 
to an unequilibrated vibration mode (Anderson-Holstein model), which can be mapped to an ef- 
fective three-state problem. Performing an iterative summation of real-time path integral (ISPI) 
expressions, we accurately reproduce known analytical results in three different limits. We then 
study the crossover regime between those limits and show that the Franck-Condon blockade persists 
in the quantum-coherent low-temperature limit, with a nonequilibrium smearing of step features in 
the IV curve. 



PACS numbers: 73.63.-b, 72.10.-d, 02.70.-c 

Introduction. — Understanding quantum transport in 
nanoscale electronic systems with vibrational or mechan- 
ical ("phonon") degrees of freedom is of topical interest 
in several areas of physics, including molecular electron- 
ics [HIl], inelastic tunneling spectroscopy |3], nanoelec- 
tromechanical systems [41, break junctions [5], and sus- 
pended semiconductor or carbon-based nanostructures 
[SHS]. The electron-phonon interaction allows to observe 
a rich variety of intriguing phenomena, such as negative 
differential conductance, the Franck-Condon blockade of 
transport, rectification, vibrational sidebands or steplike 
features in the current-voltage {IV) characteristics, and 
current-induced heating or cooling. Already the sim- 
plest nonequilibrium "Anderson-Holstein" (AH) model, 
where the nanostructure corresponds to just one spin- 
less electronic level coupled to a single oscillator mode, 
captures much of this richness |10l [TT] : for a review, see 
Ref. [3_. Analytical approaches allow to understand the 
AH model in various corners of parameter space, but 
no controlled approximation, let alone exact solution, 
connecting these corners seems in reach. One may ex- 
pect that a unified picture is available from numerics. 
However, numerical renormalization group |12| or quan- 
tum Monte Carlo (QMC) calculations [131 [H] are usu- 
ally restricted to equilibrium. For the nonequilibrium 
AH model, Han |15J employed an imaginary-time QMC 
approach followed by a double analytical continuation 
scheme; unfortunately, the latter step is plagued by in- 
stabilities |T^. A promising avenue for the AH model 
has recently been suggested by real-time path-integral 
QMC simulations |17[ [T5] , where one directly computes 
the time-dependent current. Such calculations have to 
deal with the infamous dynamical sign problem at long 
times, but in several parameter regions, especially when a 
secondary phonon bath is present, the stationary steady- 
state regime can be reached. 

In this work, we formulate and apply an alternative nu- 
merical approach, which in practice is useful unless both 
the temperature T and the bias voltage V are small. It 



is also based on a Keldysh path-integral formulation but 
does not involve stochastic sampling schemes and thus 
remains free from any sign problem. To that end, we ex- 
tend the "iterative summation of path integrals" (ISPI) 
technique jTH] to the AH model. Technical aspects of the 
present approach, in particular our mapping to an effec- 
tive three-state system via the "spin-1 Hirsch-Fye trans- 
formation" in Eq. ^ below, should also be of interest 
to QMC schemes [TJ. In essence, the ISPI method ex- 
ploits that time correlations of the auxiliary three-state 
Keldysh variable, which arise after functional integration 
over the phonon and the (dot and lead) fermion degrees 
of freedom, can be truncated beyond a certain memory 
time Tra when either T or V is finite. Together with a 
convergence scheme designed to eliminate systematic er- 
rors due to the finiteness of , such calculations allow to 
obtain numerically exact results. The ISPI method has 
already been successfully applied to the spinful Ander- 
son model |19ff2T| , where instead of the phonon a local 
charging interaction is present. While we focus on the 
simplest version of the AH model with a single unequi- 
librated [22 phonon mode here, the conceptual gener- 
alization to include Coulomb interactions, more phonon 
modes, or several dot levels is straightforward. We bench- 
mark our ISPI code against three different analytical ap- 
proaches and then study the crossover between the re- 
spective regimes. 

AH Model. — We consider the AH Hamiltonian, H = 
Hjn + Ht + Hi, describing a molecular level with tun- 
nel coupling {Ht) to metallic source and drain contacts. 
Taking a single spinless dot level (fermion annihilation 
operator d) with energy e and a boson mode (annihila- 
tion operator b) of frequency f2, the isolated molecule 
Hamiltonian is [3 (we use units with h = kg = 1) 

H,r, = nb%+[e + X{b + h^)]nd (1) 

with rid = d^d and the electron-phonon coupling strength 
A. The lead Hamiltonian Hi is taken in the standard 
wide-band approximation [23] , with relaxation processes 
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assumed fast enough to have Fermi functions with tem- 
perature T as electronic distributions; their chemical po- 
tential difference defines the bias voltage V . The tun- 
nel coupling then introduces the hybridization energy 
scales Tl and for the left/right lead |23]. For sim- 
plicity, we focus on the symmetric case in what follows, 
Vl ~ Ffl — r/2. As observable of main interest, we 
study the steady-state current / through the molecule. 

Analytical approaches. — Before turning to a descrip- 
tion of the ISPI scheme, let us briefly summarize the 
analytical approaches to the AH model that we employ 
to benchmark our method, (i) For A/F ^ 1, perturba- 
tion theory in the electron-phonon coupling applies and 
yields a closed IV expression for arbitrary values of all 
other parameters We note that the solution of the 
AH model with a very broad dot level |25l E5] corre- 
sponds to this small-A regime, (ii) For high temperatures, 
T 3> F, a description in terms of a rate equation is pos- 
sible [23] . We here use the simplest sequential tunneling 
version with golden rule rates )27| . For small A, the cor- 
responding results match those of perturbation theory, 
while in the opposite strong-coupling limit, the Franck- 
Condon blockade occurs and implies a drastic current 
suppression at low bias voltage [8l[28]. (iii) For small os- 
cillator frequency, 57 ^ min(r, eT^), the nonequilibrium 
Born-Oppenheimer (NEBO) approximation is controlled 
and allows to obtain / from a Langevin equation for the 
oscillator "W]. For small A, this approach is also con- 
sistent with perturbative theory, while for high T, NEBO 
and rate equation results are found to agree. 

Keldysh path-summation formulation. — We now start 
from the textbook time-discretized coherent-state repre- 
sentation of the Keldysh generating functional pTl 152] . 
The short-time propagator on the forward/backward 
branch of the Keldysh contour, e^*''*^, where St denotes 
the discrete time step, then allows for a Trotter breakup, 
^TiStH ^ ^^iStHi^TiSt{H-Ht) ^ ^-^^-^ the Systematic er- 
ror in observables scaling ^ It is useful to choose 
Hi = H„i — flb^b, see Eq. ([T|, where the auxiliary rela- 
tion 

allows to effectively decouple the electron-phonon inter- 
action in terms of a three-state variable = 0, ±1 de- 
fined at each (discretized) time step tj along the for- 
ward/backward (cr = ±) part of the Keldysh contour, 
where rj = {tj,a). Below, we also use the notation 
77 ± 1 = (ij±i,o') with periodic boundary conditions on 
the Keldysh contour. It is crucial for the construction of 
the coherent-state functional integral that Eq. ^ is nor- 
mal ordered. The "spin" variable Sjj picks up the three 
terms in Eq. ([2| and acts like a Hubbard-Stratonovich 
auxiliary field, similar to the Ising field employed in the 
Hirsch-Fye formulation of the Anderson model |191 155] . 
The bosonic (phonon) scalar field and the fermionic (dot 
and lead electrons) Grassmann fields appearing in the 



Keldysh path integral are then effectively noninteracting 
but couple to the time-dependent auxiliary spin variable. 
Hence those fields can be integrated out analytically and 
the time-dependent current, I{tj), follows from a path- 
summation formula for the generating functional, 

Z= deti?[{4], (3) 

{s^=0,±l} 

where the matrix Dr/n' (in time and Keldysh space) de- 
pends on the complete spin path {s}. Specifically, we 
obtain D = —iB{G^^ — S), where G^^ is the discretized 
inverse Green's function of the dot as in Refs. [1^ [ST] 
but with the modified spin-dependent matrix elements 
[~iG'^^] ^ = —s-q- The self-energy matrix S de- 
scribes the traced-out leads. We find S,,,,' ^ only when 
Sr, = ±1, where it coincides with the usual (wide-band 
limit) expression |23| . Finally, the diagonal matrix B 
(quoted here for e = 0) with 

B,„ = As^e^^''^' (4) 

encapsulates all phonon effects, where Gph is the dis- 
cretized phonon Green's function, see Ref. [31], and we 
used the notation = 1 and A±i = ±(l/2)e~'^''''*/^. 
By including a source term in D, it is straightforward 
to numerically extract the time-dependent current I{tj) 
|19] . For sufficiently long times tj, I{tj) reaches a plateau 
yielding the steady-state current of interest. 

ISPI implementation. — Starting from the formally ex- 
act path-summation formula [Eq. ([3|], the ISPI algo- 
rithm can now be adapted from Ref. |19 : For finite T 
or V , matrix entries D,,,,' involving large time differences 
|t — i'l are exponentially small. We put these to zero be- 
yond a memory time Tm = K5t, where K denotes the 
number of time slices kept in the memory. The numeri- 
cal computation of the memory-truncated path summa- 
tion in Eq. (|3| is then possible in an iterative way |T9| 
without additional approximations and, for given 5t and 
K, yields the steady-state current I{6t,K). While this 
formulation is exact for St ^ and sufficiently large 
K (long memory time), K cannot be chosen arbitrarily 
large in practice and an extrapolation scheme is neces- 
sary |19| . Convergence of the extrapolation requires suffi- 
ciently high T or V , for otherwise the necessary memory 
times are exceedingly long. For the results below, we 
used K < 4 and 0.3 < TSt < 0.35. The shown current 
follows by averaging over the i5(-window, with error bars 
indicating the mean variance. Additional ISPI runs for 
0.18 < rSt < 0.22 and 0.3 < TSt < 0.4 were consistent 
with these results, and we conclude that small error bars 
indicate that convergence has been reached. For typical 
parameters and K — 4, our ISPI code yielding I{St,K) 
runs for « 11 CPU hours on a 2.93 GHz Xeon processor. 

Benchmark checks. — Next we show that the numerical 
ISPI results are consistent with analytical theory for the 
IV curves in all three parameter limits mentioned above. 
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Figure 1. (Color online) Current / (in units of eT/h) vs bias 
voltage V (in units of F/e) for A = O.Sr, = T, e = 0, 
and T = r. The ISPI data are depicted as filled red circles, 
where the dotted red curve is a guide to the eye only and the 
error bars are explained in the main text. We also show the 
results of perturbation theory in A (dashed black curve) and 
of the rate equation (solid blue curve). The upper (lower) 
inset shows the corresponding result for T = 3T (T — F/S). 
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Figure 2. (Color online) Same as Fig. [T] but for Q — 0.5F 
(main panel) and O = 2r (inset), both for T = P. 



For clarity, we focus on a resonant level with e = here. 
Let us then start with the case of weak electron-phonon 
coupling, A = 0.5r. Figure [T] compares our ISPI data 
for = r to the respective results of perturbation the- 
ory in A and of the rate equation. As expected, for this 
parameter choice, perturbation theory essentially repro- 
duces the ISPI data. The rate equation is quite accurate 
for high temperatures, but quantitative agreement with 
ISPI was obtained only for T > lOP. Note that the ISPI 
error bars increase when lowering T due to the growing 



Figure 3. (Color online) Same as Fig.[T]but for f2 — 0.5F and 
A = F. The main panel is for T = F and compares the ISPI 
results to NEBO predictions. The insets are for T — 3T and 
T = F/3, respectively, where also the rate equation results 
are shown. Notice that in contrast to ISPI, the rate equation 
predicts an unphysical current blockade for T — F/S. 

memory time (t^) demands. The effect of changing the 
phonon frequency CI is illustrated in Fig. [2] taking T = T 
but otherwise identical parameters. Again perturbation 
theory is well reproduced. Next, Fig. [3] shows ISPI re- 
sults for a slow phonon mode, f2 = r/2, with stronger 
electron-phonon coupling, A = P. In that case, pertur- 
bation theory in A is not reliable and the rate equation 
is only accurate at the highest temperature (T = 3r) 
studied, cf. the upper left inset of Fig. [3] However, we 
observe from Fig. [3] that for such a slow phonon mode, 
NEBO provides a good approximation for all tempera- 
tures and/or voltages of interest. We conclude that the 
ISPI technique is capable of accurately describing three 
different analytically tractable parameter regimes. 

Franck- Condon (FC) blockade. — Next we address the 
limit of strong electron-phonon coupling A, where the 
rate equation approach yields a FC blockade of the cur- 
rent for low bias and T ^ P [55] . Sufficiently large A can 
be realized experimentally, and the FC blockade has in- 
deed been observed in suspended carbon nanotube quan- 
tum dots [S]. For an unequilibrated phonon mode with 
intermediate-to-large A, understanding the FC blockade 
in the quantum-coherent low-temperature regime, T < P, 
is an open theoretical problem. Here multiple phonon 
excitation and deexcitation effects imply a complicated 
(unknown) nonequilibrium phonon distribution function, 
and the one-step tunneling interpretation in terms of FC 
matrix elements between shifted oscillator parabolas [25] 
is not applicable anymore. We here study this question 
using ISPI simulations, which automatically take into ac- 
count quantum coherence effects. 

In Fig. |4] the crossover from weak to strong electron- 
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Figure 4. (Color online) ISPI data for the IV curves from 
weak (A = O.SF) to strong (A = 4r) electron-phonon coupling, 
with n = 2r. The main panel is for T = 0.2r, the inset for 
r = r. We used a dense voltage grid yielding smooth IV 
curves. Error bars are not shown but remain small, cp. Fig.[l] 



phonon coupling A is considered. The inset shows IV 
curves for T = F, where we observe a current blockade for 
low voltages once A > 2F. The blockade becomes more 
pronounced when increasing A and is lifted for voltages 
above the polaron energy A^ /fi [5S] . Remaxkably, the FC 
blockade persists and becomes even sharper as one enters 
the quantum-coherent regime (here, T = 0.2F), despite 
of the breakdown of the sequential tunneling picture. We 
also observe a nonequilibrium smearing of phonon step 
features in the IV curves in Fig. [4] cf. also Refs. |B] |2S|- 

Conclusions. — We have extended the iterative simula- 
tion of path integrals (ISPI) technique to the Anderson- 
Holstein model, which is the simplest nonequilibrium 
model for quantum dots or molecules with an intrinsic 
bosonic (phonon) mode. Our formulation exploits a map- 
ping to an effective three-state system and reproduces 
three analytical theories valid in different parameter re- 
gions. This extension of the ISPI approach then captures 
the full crossover between those limits unless both T and 
V are very small. For strong electron-phonon coupling 
and an unequilibrated phonon mode, we find that the 
Franck-Condon blockade becomes even more pronounced 
as one enters the quantum coherent regime. 
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